arXiv:1505.00793v2 [cond-mat.soft] 6 May 2015 


Crystal nuclei in melts: A Monte Carlo simulation of a model for attractive colloids 


Antonia Statt 1,2 , Peter Virnau 1 and Kurt Binder 1 
1 Institut fur Physik, Johannes Gutenberg- Universitat Mainz, Staudinger Weg 9, D-55099 Mainz, Germany 
2 Graduate School of Excellence Materials Science in Mainz, Staudinger Weg 9, 55128 Mainz, Germany 

As a model for a suspension of hard-sphere like colloidal particles where small nonadsorbing 
dissolved polymers create a depletion attraction, we introduce an effective colloid-colloid potential 
closely related to the Asakura-Oosawa model but that does not have any discontinuities. In simu¬ 
lations, this model straightforwardly allows the calculation of the pressure from the Virial formula, 
and the phase transition in the bulk from the liquid to crystalline solid can be accurately located 
from a study where a stable coexistence of a crystalline slab with a surrounding liquid phase oc¬ 
curs. For this model, crystalline nuclei surrounded by fluid are studied both by identifying the 
crystal-fluid interface on the particle level (using suitable bond orientational order parameters to 
distinguish the phases) and by “thermodynamic” means. I.e., the latter method amounts to com¬ 
pute the enhancement of chemical potential and pressure relative to their coexistence values. We 
show that the chemical potential can be obtained from simulating thick films, where one wall with 
a rather long range repulsion is present, since near this wall the Widom particle insertion method 
works, exploiting the fact that the chemical potential in the system is homogeneous. Finally, the 
surface excess free energy of the nucleus is obtained, for a wide range of nuclei volumes. From this 
method, it is established that classical nucleation theory works, showing that for the present model 
the anisotropy of the interface excess free energy of crystals and their resulting nonspherical shape 
has only a very small effect on the barrier. 


INTRODUCTION 

Understanding how crystalline solids form from fluid 
phases is a problem that both poses significant intellec¬ 
tual challenges and is a basic input for materials science 
m- The standard picture implies that nanoscopic nu¬ 
clei of the solid are formed due to spontaneous thermal 
fluctuations when the fluid density has a value exceeding 
the fluid density at fluid-solid coexistence. Due to the 
cost in surface excess free energy caused by the crystal- 
fluid interface HE), a free energy barrier in phase space 
needs to be overcome in every event of this so-called ho¬ 
mogeneous nucleation [DE1I1I3 process. In the present 
paper, we are neither concerned with the further crystal 
growth processes m o uni that follow after the crystal 
nuclei have been formed, nor shall we consider the fact 
that in practice heterogeneous nucleation 0EEIHH] may 
dominate (since the free energy barrier may be signif¬ 
icantly reduced for nuclei attached to the walls of the 
container, or to seed particles E EO EH ED E3 , etc.). 

Due to the large value of the nucleation barrier (which 
in typical cases is several ten’s of the thermal energy 
knT | 1 , [H ':5J [7]) and the nanoscopic size of these nu¬ 
clei (typically they contain only a few hundred particles 
EJ H El 13) nucleation events are rare and difficult to 
observe directly (either in experiment or in simulation). 
Also the detailed properties of the nuclei that are formed 
are non-universal and depend upon the particular sys¬ 
tem under study. Experiments often infer the nucleation 
rate indirectly from later stages of the process, when the 
phase transformation from fluid to solid reaches comple¬ 
tion [SMI. Of course, this procedure inevitably involves 
questionable assumptions. Simulations of nucleation ki¬ 
netics, however, often are feasible only under conditions 


where the barrier is only of order 10k bT or less E1E3- 
Then it is inevitable that the nuclei are extremely small 
and strongly fluctuating objects, for which the “classical” 
description in terms of a competition between bulk and 
surface free energies [3IHEl 13 fails [HEE]. 

Many of the difficulties alluded to above are already en¬ 
countered when one considers nucleation of fluid droplets 
from a supersaturated vapor mm usi- However, for nu¬ 
cleation of crystals from melts there occur three addi¬ 
tional complications (i) the interface tension (^y (fi)) be¬ 
tween a crystal surface and the fluid depends on the ori¬ 
entation of the unit vector n normal to the surface rela¬ 
tive to the axes of the crystal lattice Emn]. The shape of 
the nucleus then is not a spherical droplet as in the case 
of vapor to liquid nucleation, but also exhibits anisotropy. 
For strong enough anisotropy, facetted nanocrystals are 
expected EHHMI- I n general, the equilibrium crystal 
shape is nontrivial and can only be obtained from 7 (n) 
via the Wulff construction pUH55] . (ii) Also kinetic pro¬ 
cesses, such as attachment of particles to the surface of a 
nanocrystal, may depend on the surface orientation [lO]. 
If the nucleation barrier is not very high, and nucleation 
plus growth proceed relatively fast, crystal shapes rather 
different from the equilibrium shape could result 0 . (hi) 
Often, a fluid at the melting/crystallization transition is 
rather dense, and enhancing the density further a glass 
transition is encountered [3]. Then the variation of the 
kinetic prefactor in the nucleation rate with density be¬ 
comes an important factor 0E3- Far from the equi¬ 
librium coexistence conditions a separation of a static 
Boltzmann factor involving the nucleation barrier and a 
kinetic prefactor may be impossible. Similar problems 
arise in other systems where the liquid to solid transition 
is not driven by enhancing the density, as is done in the 
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generic case of colloids described as a hard sphere fluid 
@123 EMI, but by reducing the temperature: in flu¬ 
ids such as silica @ mi B 2 I, triphenyl phosphite 1431 l44j 
and ortho-terphenyl strong supercooling is possible 
and ultimately a glass transition Hi 2 is reached. Nu- 
cleation competing with glass formation [3] shall not be 
considered here, and also the interplay of crystallization 
with a liquid-liquid unmixing, as it has been proposed 
for water and various other liquids [HI 0 B] is out of our 
focus. 

The present paper hence is only concerned with the 
first problem, of predicting the nucleation barrier associ¬ 
ated with the formation of crystalline nuclei, but avoiding 
the assumption that 7 (ft) actually is isotropic and the nu¬ 
clei are spherical. A central element of the present paper, 
however, is, that the nucleation barrier can be obtained 
directly from a careful analysis of the equilibrium con¬ 
ditions for the nucleus and the surrounding fluid phase 
S3- While in the thermodynamic limit this surround¬ 
ing fluid phase is metastable, we show that in finite sys¬ 
tems there actually occurs a stable equilibrium between 
the small crystalline nucleus and this surrounding fluid 
phase. We show, using a simple model for a suspension 
of colloids with short-range attractive interactions, that 
in this method there is no longer any need to obtain 7 ( 71 ) 
and to carry out the Wulff construction. While the pre¬ 
cise shape of the crystalline nucleus is not found in this 
approach, nucleation barriers can be estimated with rea¬ 
sonable precision. 

In the next section, we shall describe the model that we 
use, a variant of the Asakura-Oosawa (AO) model [451 - 
[50] for colloid-polymer mixtures EDE3, and summarize 
what is known on the thermodynamic properties of this 
model in the bulk. We recall here that the AO model 
contains the widely studied case of the hard sphere fluid 
as special case, but has the merit that the control pa¬ 
rameters of the attractive part of the potential (its depth 
and its range, respectively) can easily be tuned (experi¬ 
mentally this could be done by varying the concentration 
of polymers in the suspension, and their size via choos¬ 
ing different molecular weights ED E2])- Since in this 
way both the density gap p m — Pf (onset of melting of 
the crystal occurs at a density p = p m , onset of freezing 
of the fluid occurs at p = pf) and the interfacial tension 
7 (n) can be varied over a wide range, the system qualifies 
as a model suitable to stringently test the theory. 

In section 3, the preparation of the inhomogeneous 
systems, where a crystalline nucleus coexists with sur¬ 
rounding fluid in a finite simulation box, is described in 
detail, since the equilibration of such systems exhibiting 
two-phase coexistence is notoriously difficult EDIESES- 
[58] . We also introduce the local bond orientational order 
parameters m [60] that are used to identify (in each 
system configuration that needs to be further analyzed) 
where the crystalline nucleus (and its interface separat¬ 
ing it from the surrounding fluid) is located. Note that 
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FIG. 1: Potential ?7(r)/fcflT versus distance r between two 
colloidal particles, for the SoftEffAO model, and compared to 
U £ ff( r )/U bT for the EffAO model. In both cases, q = 0.15 
and r/p = 0.1 is chosen, and a c = 1 is used as unit of length. 

the identification of interfaces between coexisting phases 
at the single particle level involves important ambiguities 
RHH5S] ; we feel that previous work where the estimation 
of the nucleation barrier as function of the nucleus size 
was attempted, based on such methods, possibly some¬ 
what suffered from this ambiguity. 

In Sec. 4 the thermodynamic method utilized in the 
present paper, where the density pressure p^, and 
chemical potential pg of the liquid surrounding the crys¬ 
tal nucleus are accurately estimated, for a wide range of 
total particle densities p (varying the simulation box vol¬ 
ume Vbox for several choices of the number N of colloidal 
particles). We explain how from an analysis of this equi¬ 
librium we can infer both the volume V* of the crystalline 
nucleus, and the associated free energy barrier A F*. 

Finally, in Sec. 5 we give a brief summary of our re¬ 
sults, and an outlook on further related problems. 

THE MODEL AND ITS BULK PROPERTIES 

From the Asakura Oosawa (AO) model to the Soft 
Effective Asakura-Oosawa (Soft EffAO) model 

Colloid-polymer mixtures [50)452] contain three ingre¬ 
dients: colloidal particles, hard sphere-like in a size range 
from 10 nm to 1 pm; flexible polymers, which have a high 
enough molecular weight to form coils in solution with a 
radius in the range from 10 nm to lOOnm; and a small 
molecule fluid in which the colloidal particles are sus¬ 
pended, and which must be chosen such that it acts as 
a good solvent for the polymers or a Theta solvent, in 
which the polymers take a random walk-like structure 
m- By suitable coating of the colloids with a surfactant 
layer one can avoid polymer adsorption on the colloids, 
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and screen out the attractive van der Waals interaction 
between the colloids, avoiding their aggregation. 

The AO model provides a simplified description of the 
statistical thermodynamics of such a system, represent¬ 
ing the colloids as hard spheres of diameter a c , and the 
polymers as soft spheres of diameter a p . While polymer- 
colloid overlap is also strictly forbidden, polymers can 
overlap each other with zero energy cost FlMoUj . As is 
well known the depletion effect causes an (entropic) ef¬ 
fective interaction between the colloids. Note that the 
molecular solvent in this model is not explicitly regarded 
at all. If the size ratio q = <r p /cr c is small enough 
(q < q* = 2/>/3 — 1 = 0.1547 [SS]), one can derive an 
equivalent one-component model (which we refer to as 
“Effective AO model” in the following ) of a colloidal 
fluid, where the particles interact with the following ef¬ 
fective potential: 

U e ff(r < cr c ) = oo, 

/3U e ff(a c < r < a c + a p ) = ^ a p z p (l + g -1 ) 3 - (1) 

( _ 3r/a c (r/ffc ) 3 \ 

V 2(1 + 5) 2(1 + qf) ’ 

and U e ff(r > a c + o>) = 0. (2) 


Here r is the radial distance between two colloidal parti¬ 
cles, and z p = ey^p^p/ksT) /A 3 the polymer fugacity (/i p 
is the polymer chemical potential and A p their thermal 
de Broglie wavelength). In the following we shall denote 
the model defined by Eq. [l] where the polymer degrees 
of freedom were integrated out and hence all informa¬ 
tion on the distribution of polymers is lost, as Effective 
Asakura-Oosawa (EffAO) model. 

While Eq. [l] is convenient for analytical work (e.g. us¬ 
ing density functional theory [69]), it is cumbersome to 
use in simulations: the singularity of U e ff(r) at r = a c 
precludes the use of the Virial formula ED] to compute 
the pressure in the system (alternative methods to com¬ 
pute the pressure in systems containing hard particles ex¬ 
ist ogi but are difficult to apply and require huge nu¬ 
merical effort; also the use of Molecular Dynamics meth¬ 
ods mm would be very inconvenient). 

Since real colloidal particles anyway never behave as 
strictly hard spheres EH Eg, however; there is no need 
to stick exactly to the potential, Eq. |T] but one can 
equally well consider a model where the singular behav¬ 
ior at a = r c has been replaced by a smooth function. 
Thus we replace the hard core part of Eq. [T] by the fol¬ 
lowing model, called the Soft Effective Asakura-Oosawa 
(SoftEffAO) model, while the other parts of Eq. [l] are 
maintained, 


Urep{p ) — 
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The constants b, e are chosen as b = 0.01 and e = 
0.98857; of course, the form chosen for Eq. [3] is rather ar¬ 
bitrary, other functions could be used equally well. Figjl] 
gives a comparison of both potentials, Eq. [I] and Eq. [3] 



FIG. 2: Normalized pressure pu^/UbT plotted vs. T) for 
both the EffAO model (squares) and the soft EffAO model 
(circles). Full curves through the data points are guides to 
the eye only, and statistical errors are smaller than the size of 
the symbols, and hence not shown. Note that the equation of 
state contains two disjunct functions; one for the liquid branch 
(pe{v)) and one for the crystalline solid (p c (p)) (which has a 
face-centered cubic (fee) structure), respectively. Note that 
the data for the liquid branch comprise also metastable states 
(for p > Pcoex ) and the data for the solid branch comprise 
metastable states as well (for p < p CO ex)- The coexistence 
pressure p CO ex is found from the interface velocity method 
(see Fig. 3), and the coexistence packing fractions then result 
from r/f = r/e{p C oex),ri,Ti = r/c(Pcoex), with rie{p),Vc(p) being 
the inverse functions of pe(ri),Pc(v)- The result for the liquid 
were obtained using N = 4000 colloidal particles and using 
both the NpT ensemble and the NVT ensemble (see text) to 
show that finite size effects are negligible. 


For this SoftEffAO model the application of the Virial 
formula is straightforwardly possible, and as shown in 
Fig- [2], the equation of state of both models is very sim¬ 
ilar, for the choice q = 0.15, and r] p = 0.1 that will be 
studied here. Here we have defined the packing fraction p 
of the colloids and the polymer reservoir packing fraction 
T] r p as (p = N/V box ) 

V = (kct 3 c /6)p, rf p = (tt<j 3 /6) exp (p p /k B T). (4) 

Note that p” is defined such that it would give the 
packing fraction of polymers (?y” = {Tra p /6)N p /Vbo X ) if 
the simulation box would contain the ideal gas of poly¬ 
mers only. 


Bulk phase behavior of the SoftEffAO model 

The pressure isotherm shown in Fig. [2] contains two 
branches, the liquid branch (at the left) and the solid 
branch (representing a face-centered cubic (fee) crystal, 
at the right). Actually the solid branch was computed in 
the NpT ensemble, choosing a cubic box with N = 4000 
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particles, arranged at the sites of a fee lattice as ini¬ 
tial condition. Simulations where only the total volume 
V = L x L y L z is allowed to fluctuate or where the simula¬ 
tion box linear dimensions L x ,L y ,L z in the three carte¬ 
sian directions are allowed to fluctuate independently 
gave (within statistical errors) identical results for the 
function 77 = 77 ( 77 ) that is shown in Fig. [2] In our sim¬ 
ulations, we have disregarded the alternative hexagonal 
close-packed (hep) structure throughout, since it is not 
expected to be the stable phase [5DH55] . Of course, it 
nevertheless would be interesting to study crystal nuclei 
having the hep structure surrounded by fluid, applying 
the technique of the present work, and to compare their 
nucleation barriers to those of fee nuclei. However, this 
interesting issue must be left to future work. 

The liquid branch was obtained from both runs in the 
NpT ensemble and in the NVT ensemble (sampling then 
p at given fixed 77 from the virial expression [7(71 l73l ITT ] ). 
Within statistical errors, the results from the NpT and 
NVT ensembles agree precisely, indicating hence that fi¬ 
nite size effects on the isotherms for N = 4000 are already 
negligible. We also note that for the ranges of 77 shown 
we never could observe a spontaneous crystallization of 
the liquid, nor a spontaneous melting of the solid. This 
implies that the time span of the simulation runs was 
too short to observe a decay of a metastable phase via 
nucleation and growth. 

It also is interesting to ask for the effect of choosing 
different parameters b and e in the potential, Eq. [3] One 
can show that increasing b has the effect of displacing the 
liquid branch slightly (and almost uniformly) to the left. 
The effect on the solid branch is much more pronounced, 
however, it is shifted to the left by much larger extent, 
and hence the width of the liquid-solid coexistence region 

decreases 1223- 

In thermal equilibrium, the liquid phase and the crystal 
in Fig. [2] can coexist at a single pressure only, p = p CO exi 
where also the chemical potentials pe(p), p, c (p) of both 
phases are equal, 


/fi (Pcocx ) — Pc.iPc.oex ) — Me 


(5) 


But from Fig. [5] we see that there occurs an extended 
region of hysteresis, where either the liquid ( p > p coex ) 
or the crystal (p < p CO ex ) is metastable in our simula¬ 
tion. Hence the accurate estimation of p coe x is a nontriv¬ 
ial problem. We follow the procedure of Zykova-Timan 
et al. [65] by locating p coex from studies of “slab config¬ 
urations” (Fig. [ 3 ]) of liquid-solid coexistence in the NpT 
ensemble. We estimate the interface velocity from the 
time derivative of the volume ( dV/dt) as a function of 
pressure p , identifying p coex from the pressure for which 
(dV/dt) changes sign (Figs. EH- The initial system 
at each pressure is chosen as rectangular parallelepiped, 
V = L x L x L z , choosing L z « 5L, and a crystal slab of 
extension L c (« L z / 2) in the z-direction, with two (001) 
interfaces (which for the length of the runs performed 



FIG. 3: Snapshot picture of a EffAO configuration of a sys¬ 
tem with linear dimensions L x L x L z with L z = 5 L, the 
crystal slab on the left side and the liquid on the right side. 
Here n = 6 and N = 4000 was chosen, appropriate for p = 8.0. 
Particles belonging to the liquid are shown in blue color, par¬ 
ticles belonging to the crystal are shown in red color, while 
particles in the center of the interface are displayed in green 
color. For the cutoff values for qe and <74 refer to [82) . 


are still non-interacting, even if the crystal gradually 
melts (and therefore the volume change AV increases 
with time, Fig. [4]). At each p , the linear dimension L is 
chosen as L = nd(p), n being an integer (5 < n < 12) 
and d(p) the equilibrium lattice constant at the pressure p 
that we wish to test; d(p) can be estimated from the crys¬ 
tal branch in the p vs. 77 isotherm (Fig. [3|. Also L c was 
chosen commensurate with an integer number of crystal 
planes. Such a choice of linear dimensions commensu¬ 
rate with the crystal structure is essential, since we use 
periodic boundary conditions throughout, and one must 
avoid strain deformations of the crystalline slab. The re¬ 
maining volume Ly.Lv, ( L z — L c ) contains liquid with 
a density chosen according to the liquid branch on the 
isotherm 77 ( 77 ) for the pressure that is tested (the precise 
value of L z for the initial configuration is chosen such 
that the particle number N is an integer, of course). To 
prepare the initial configuration for the NpT simulation 
of this slab configuration we actually carry out first a run 
in the NVT ensemble with V fixed through the initial val¬ 
ues of L and L z , to achieve that the liquid configurations 
in the crystal-liquid interfacial region also equilibrates 
(only particles in the liquid region take part in this ini¬ 
tial equilibration, carrying out about 10 10 Monte Carlo 
moves, i.e. displacements of single particles). 

In the simulations in the NpT ensemble then all par¬ 
ticles are considered for displacement moves, and in the 
volume moves all three linear dimensions L x ,L y and L z 
are allowed to fluctuate independently (unlike [65] where 
L x = L y = L = nd{p) was held fixed, so only L z could 
fluctuate), but keeping n fixed. This procedure required 
an order of magnitude more statistical effort than used 
in [55]. The length of the runs was 3 x 10 4 MC cycles 
(one cycle consisted of N particle displacement moves 
and one volume move), and in each case 100 indepen¬ 
dent runs were made. 

Fig. [3] shows that AV increases with time for low pres¬ 
sures (as expected, due to the progressing melting of the 
crystalline slab the volume fraction of the liquid that re¬ 
quires more volume increases) and for high pressures the 
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FIG. 4: Volume change AV as a function of Monte Carlo 
(MC) cycles for different normalized pressures pa^/kgT, as 
indicated in the figure, for the effective AO model (a) and 
for the soft effective AO model (b). The system contains 
N = 9250 colloids in each case, corresponding to n = 8 lattice 
planes in x and y direction. The solid straight lines are fits, 
from which (A V/dt) is determined. Each curve is averaged 
over 100 independent runs. 


opposite trend is observed. Clearly, these data still are 
affected by statistical fluctuations, and there is also sig¬ 
nificant scatter in the variation of the interface velocities 
with pressure (Fig. [5}c). These difficulties are a serious 
limitation for the accuracy with which p coex can be deter¬ 
mined. Taking also the trend with increasing number n 
of planes (and hence increasing linear dimension L) into 
account (Fig. [5 ]d), we conclude 

p coex = 8.09 ± 0.06 (Eff AO model), 

Pcoex = 8-45 i 0.04 (soft Eff AO model). (6) 



l/n 2 

FIG. 5: (a) Velocities of volume changes (A V/dt), as deter¬ 

mined in Fig. 4, plotted vs. p, for different choices of n, as 
indicated. Results for the effective AO model and for the soft 
effective AO model (b) are shown. Straight line fits of the data 
are shown, from which the transition point ((A V/dt) = 0, 
highlighted by a horizontal straight line) p = p coex is esti¬ 
mated. (c) Pressure p where (A V/dt) = 0 plotted vs. l/(n 2 ), 
where n is the number of lattice planes in x and y direc¬ 
tion in the box. Disregarding the smallest size, the remaining 
data are fitted empirically to straight lines, yielding the result 
Pcoex = 8.09±0.06 for the EffAO model andp CO em = 8.45±0.04 
for the soft EffAO model. 


In our preliminary publication m, slightly different 
values were quoted, p coe x = 8.06 ± 0.06 ( EffAO model) 
and Pcoex = 8.44 ± 0.04 (soft Eff AO) relying on the em¬ 
pirical observation that an extrapolation against 1 /vN 
seemed to work best m- Note, however, that there is 
no theoretical justification for an extrapolation against 
i /Vn, one might rather expect an extrapolation against 
the inverse interface area should be used, as done in Fig. 
5c. As a caveat, we note that if we would have chosen 
only the tree smallest choices of n, one would have er¬ 
roneously concluded that finite size effects on pcoex are 
very small, resulting in estimates significantly off. As 
we shall see later, this error (which is less than 1% ) 
is the most critical limitation to the accuracy that we 
can reach for the prediction of nucleation barriers; how¬ 
ever, it is by no means straightforward to obtain p coex 
with significantly smaller errors. Knowing the coexis¬ 
tence pressures the corresponding packing fractions of 
the coexisting phase are readily found from Fig. [2] For 
the soft effective AO model, this implies 

Vf = 0.495 (1) , r lm = 0.636(1) (7) 


Computing the chemical potential for dense fluids 

The chemical potential /i of the colloids is a basic quan¬ 
tity for the study of phase coexistence, since (unlike the 
pressure tensor) p is strictly spatially constant in a sys¬ 
tem in thermal equilibrium, even if the system is inhomo¬ 
geneous due to the presence of interfaces, external walls, 
etc. Thus when we consider the coexistence of a crys¬ 
talline nucleus of finite size with surrounding liquid in a 
simulation box, we have 

A* = He(pe) = Aic(Pc) > Pcoex (8) 

instead of Eq. [5j which applies in the thermody¬ 
namic limit. We expect p c > pi > p CO ex, due to the 
Laplace pressure effect. Quantification of these differ¬ 
ences /J - Pcoex, Pi - Pcoex, Pc - Pcoex as function of 
nucleus size is one of the main aims of the present paper, 
which will allow us also to estimate the nucleation bar¬ 
rier, without identification of the crystal-liquid interface 
on the particle level (Sec. 4).Therefore precise estimation 
of the chemical potential is a matter of concern. 
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FIG. 6 : (a) Illustration of the method to compute the chem¬ 
ical potential of a dense colloidal fluid, using a L x L x L z 
geometry with two walls at z = 0 and at z = L z , and pe¬ 
riodic boundary conditions in x and y directions, for the 
case L = 7, L z = 30 and four choices of particle number 
(N = 750,950,1100 and 1250, respectively, at z = 15 from 
bottom to top respectively). At z = 0 a soft repulsive poten¬ 
tial acts (Eq.[4| while at 2 = L z a hard wall acts (which causes 
significant layering). The insert shows the regions where in 
each case the chemical potential p, is computed (and demon¬ 
strates that /r is independent of z in this region in each case, as 
it should be). The bulk liquid packing fraction corresponding 
to fj, can be read off from an average over the region where the 
packing fraction profile 77 ( 2 ) is constant, i.e. for 10 < z < 20 
in this example, (b) Soft wall potential V 2 (z) given in Eq|10| 
as function of distance z in colloid diameter to the wall. 



Now it is well known that the standard Widom parti¬ 
cle insertion method [78] becomes inapplicable at liquid 
densities near the melting transitions, since the accep¬ 
tance probability for virtual particle insertions becomes 
extremely small. It was early recognized, however, that 
in principle one could circumvent this problem studying 
a system with a soft repulsive wall, since near this wall 
the density is reduced and hence the acceptance proba¬ 
bility is higher m • Since the chemical potential near the 
wall must be the same as far away from the wall (where 
the liquid density may be so high that particle insertions 
fail), we know the chemical potential also in the latter 
region. 

Of course, the accuracy of this approach is a nontriv¬ 
ial problem; the spatial extent of the region where the 
chemical potential can be “measured” covers only a small 
fraction of the system when the system is large. When 
the system is rather small, as used in ESI, the density 
shows pronounced oscillations throughout the simulation 
box (“layering”), and then it is not at all clear to which 
bulk conditions the considered system would correspond. 
As a matter of fact, to our knowledge the simple idea of 
Ref. 79 has never been convincingly validated with re¬ 
spect to its practical usefulness yet. 

Fig. [ 6 ] now compellingly shows that the method indeed 
works well, for our model. At the strongly repulsive wall 
at the right side in Fig. [ 6 ] , for z = L z , we use a Weeks- 


Chandler-Anderson-type potential 


( 0+1 \ 12 

( &wl \ 

V*-) 

v^r) 


z! = L z — z, 0 < z' < tr u ,i2 1 ^ 6 , 


(9) 


with e w i = 1 and a w i = cr c /2, and V\ (z) = 0 for z' > 
(J u ,l 2 1/,6 . At the weakly repulsive wall at the left side in 
Fig. 6 , we choose a potential of the form 


V 2 (z) = 
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(if) 


with e W 2 = 0.8, e' w2 = 0.1, a w 2 = 1.5, h = 10 -8 and 
V 2 (z) = 0 for z > z c = 2.2a W 2- Note that neither V\{z) 
nor V 2 (z) exhibit a discontinuity at their cutoff values, 
also the force is continuous there. Due to the large range 
of a W 2 and the rather smooth variation of V\ (z) the lay¬ 
ering of the packing fraction r](z) quickly reaches a hor¬ 
izontal plateau. Due to the depletion of particles near 
the soft wall (and the excess density that occurs near the 
hard wall, where a pronounced layering is observed EQl) 
the packing fraction 77 ^ that is observed in the flat re¬ 
gion of V 2 (z) differs somewhat from the average packing 
fraction fj in the simulation box. Note that 


L z 

fj = L~ 1 jr,(z)dz = N/(L 2 L z ) (12) 

0 

In the example shown in Fig. [ 6 ] , we have chosen 
fj « 0.2671, 0.3384, 0.3918 and 0.4452 (from bottom 
to top), while the corresponding estimates for rji, were 
r]b = 0.2933, 0.3690, 0.4241 and 0.4777, respectively. The 
chemical potential /i (shown in the insert) needs to be 
associated with rjb , of course, since fj differs from rib by 
wall excess corrections, 

1 1 

V = Vb+ ~j—r)w 1 + 7 - 17^2 , (13) 


with 


L z /2 


f] w \ = J ( r)(z ) - ilb)dz , ri w2 = J (■ T](z ) - rjb)dz. 


L z /2 


(14) 

The data given would allow precise estimation of these 
excess densities of the walls, 77^,1 and rj W 2 , as well, but this 
is out of the scope of the present paper. We rather em¬ 
phasize that near the left wall, where 77 ( 21 ) gradually rises 
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FIG. 7: (a) Chemical potential p (in units of fcsT) plotted 

vs. packing fraction r /, for different choices of L and L z , as 
indicated. The data labeled as “Widom” (available for not so 
large packing fractions) are obtained by the standard particle 
insertion method for homogeneous bulk systems. Arrows on 
abscissa and ordinate indicate 77 / and /. i CO ex , respectively, (b) 
Chemical potential p plotted vs. pressure p, using the same 
data as in part (a). Arrows at the abscissa and ordinate 
indicate p coe x and p coe x, respectively. 


transition point 77 = 77 / = 0.495 (|T|) from Figs. ( j2|5[ o), we 
can read off Hcoex = 4.74 ([ 5 ]) as well (Fig{7|i,b). Note that 
Fig. & shows that the relation pi(p) versus p for the liq¬ 
uid branch near p coex is very accurately represented by a 
straight line, hence validating a simple Taylor expansion 


7 r 1 

Wipe) ~ Vcoex + W- {PC-Pcoex) (15) 

6 77/ 

Here we have given the pressure an index t to empha¬ 
size that Fig. [TJd and Eq. [15] describe the pressure of the 
liquid phase and not the crystal (when we consider a crys¬ 
tal nucleus surrounded by the liquid phase, the pressure 
p c in the crystal exceeds pi). 


HOMOGENEOUS SYSTEMS VERSUS 
TWO-PHASE COEXISTENCE 

Local Identification of Phases in Terms of Bond 
Orientational Order Parameters 


to 77 b, we can use a significantly broad slab (of a width 
Az of the order A z = 1 ) where rj(z) is small enough 
that the Widom particle insertion method allows an ac¬ 
curate sampling of the chemical potential. Although 77 ( 2 ) 
is rapidly varying in these intervals of z where the sam¬ 
pling of p is done ,we do find that p is nicely constant in 
each case, the statistical fluctuations being smaller than 
the thickness of the horizontal straight lines shown in 
the insert.Although this method relies only on a small 
fraction (A z/L z ) of the volume that is simulated, never¬ 
theless a satisfying accuracy of the function p(r)b) in the 
fluid phase is reached. One could have chosen a soft wall 
potential at both walls, and thus improve the statistics 
for the estimation of the chemical potential somewhat. 
Such a choice would in particular be necessary in cases 
where complete wetting of the fluid at melting conditions 
at the hard walls occur, since then there would be an ex¬ 
tended range of strong layering. In our case, the layering 
extends over a few particle diameters only, and thus does 
not matter. 

Fig. [7] presents then a plot of p vs. rj over a wide range 
(0.3 < 77 < 0.52) of the fluid branch of the isotherm, 
including hence also its metastable part. For not too 
large values of 77(77 < 0.46) also estimates recorded for 
bulk L x L x L systems (with no walls and periodic 
boundary conditions throughout) are included, apply¬ 
ing the standard particle insertion method. These es¬ 
timates are in excellent agreement with the results found 
from the method described above, only for the last point 
(77 ~ 0.455) the result of this direct method is slightly 
off, due to a too small acceptance rate of the particle in¬ 
sertions. We also include many different choices of both 
L and L z in Fig. [7] and demonstrate that our results are 
not hampered by finite size effects. Since we know the 


When we want to analyze the phase coexistence be¬ 
tween a crystalline nucleus and surrounding liquid phase 
(the associated chemical potential then can be inferred 
from Fig. S>) as a function of the nucleus volume, it is 
mandatory to identify where inside the simulation box 
the crystalline nucleus actually is located, for each con¬ 
figuration of the colloidal particles that is analyzed. Due 
to statistical fluctuations fluid particles in the interfacial 
region may get attached to the crystal, as well as the 
inverse process. While in equilibrium there should be 
just fluctuations of the size and shape of the crystalline 
nucleus around their proper average values, these fluc¬ 
tuations will also lead to random displacements of the 
center of mass of the crystalline nucleus. In a long sim¬ 
ulation run, we hence should have a diffusive motion of 
the crystalline nucleus in our system. 

Of course, reliable “measurements” of pc can only be 
taken in such subboxes of the system which are not af¬ 
fected by the presence of the nucleus; i.e. remote also 
from the region of the crystal-liquid interface in the an¬ 
alyzed configuration of the system. Note that also the 
orientations of the lattice planes at the crystal nucleus 
can be arbitrary, they should not be biased by the axes 
of the simulation box, if the simulated system is large 
enough. 

The identification of whether particles are part of the 
liquid phase or part of the crystalline nucleus is based on 
the well-known Steinhardt bond orientational order pa¬ 
rameters [521 EDI. They are defined in terms of spherical 
harmonics as 

1 N(i) 

'y 1 Yem{0( r ij)i < P( r ij)) ( 16 ) 










FIG. 8 : (a) Distribution P(<j 4 , qe) in the liquid ( at 77 = 0.485 
), hep and fee phases (at rj = 0.65) (b) Variation of ( 94 ) and 
and (qe) with 77 in the liquid and fee phases, for the range 
from 77 = 0.53 to 77 = 0.565. (c) Distribution recorded for 
the coexistence of a crystal nucleus surrounded by fluid. This 
case refers to an average packing fraction 77 = 0.528 in the 
simulation box and 777 = 0.512 in its liquid part. 


where fij is the distance vector from the considered col¬ 
loidal particle, labeled by index i, to one of its N ( i ) neigh¬ 
bors within a cutoff radius (defining what the “nearest 
neighbors” are). The polar angles 9(fij) and <j>(fij) of 
this distance vector are chosen with respect to an arbi¬ 
trary but fixed reference frame. We use a cutoff radius 
of 1.3 <7 here, determined by the first minimum of the 
radial distribution function g(r) measured in a bulk fee 
crystal. It is advantageous [60] to average the qe m {i) over 
a particle and all its nearest neighbors 


1 JV(i) 

Qem(i) = N(i) 1 (17) 


where ki = 0 means that one takes from 

The order parameter that finally is used is then 
as 


Eq. 16 


defined 


&(*) = ^ 

' v ' m——i 


1/2 


(18) 


In order to distinguish fee crystals from liquids it is use¬ 
ful to consider the resulting distribution of <74 and qe 
(Fig. [ 8 ^). One sees that in the pure phases the dis¬ 
tribution are well separated from each other. Since it 
is conceivable that in simulations of phase coexistence 
also hexagonal closed packed (hep) regions compete with 
fee regions, also data for P(g 4 , q 6 ) for a hep crystal are 
included: However, in our study we have not encoun¬ 
tered any evidence that this structure plays a role for 


our model. It is also gratifying to note that the peak po¬ 
sition of qe and <74 are basically independent of 77 , when 
one studies fee crystals in the relevant region of packing 
fractions (Fig. |8jb). 

The situation is somewhat different when one consid¬ 
ers phase coexistence between a crystalline nucleus and 
surrounding fluid (Fig. Now one finds a “ridge” 

connecting the sharp peak (representing the crystalline 
nucleus) with the broad peak representing the liquid. 
Clearly there is an ambiguity whether particles yielding 
a “signal” in this intermediate region are to be “counted” 
as part of the liquid or as part of the crystal. We shall 
return to this problem in the next section. 


Preparation and local Characterization of Systems 
exhibiting Phase Coexistence Between a Crystalline 
Nucleus and Surrounding Fluid 

We study systems of cubic shape (L x L x L boxes 
with periodic boundary conditions) at constant volume 
and a particle number N chosen such that the pack¬ 
ing fraction 77 is inside the two-phase coexistence region, 
namely 0.525 < 77 < 0.530, for a total number of col¬ 
loids N = 6000, 8000 and 10000 and start with a situa¬ 
tion where a compact crystal nucleus (of about the right 
volume V c and packing fraction q m , so that the remain¬ 
ing particles are compatible with a surrounding liquid 
of the expected density) is initially already present. In 
the quoted parameter range, one finds that after a (very 
long) equilibration stage the system settles down at the 
stable equilibrium situation, where the crystal has some¬ 
what grown or shrunk in size and adjusted its shape, such 
that the proper equilibrium is established. 

When we try to do this at significantly smaller 77 
(but still in the two-phase coexistence region, e.g. for 
77 = 0.520), we regularly find that the crystalline nu¬ 
cleus melts, and the system ultimately becomes a ho¬ 
mogeneous fluid; when we try to do this at significantly 
larger 77 , e.g. for 77 = 0.535, we find that the compact 
droplet undergoes ultimately always a fluctuation to an 
elongated, more or less cylindrical, shape, connected into 
itself through the periodic boundary condition. For still 
larger 77 (0.54 < 77 < 0.57), even the cylinder is unsta¬ 
ble, and the system develops towards a “slab” configu¬ 
ration, where a crystal domain extending in two direc¬ 
tions throughout the box is separated from the liquid 
by two planar interfaces (in this region one must work 
with L x L x L z boxes, varying 77 by changing L z and 
choosing L commensurate with an integer number of lat¬ 
tice planes, of course, to avoid that the crystal domain is 
strained, cf. Sec. 2.2). Thus, there is only a rather small 
range of packing fractions where the desired equilibrium 
between a crystalline nucleus and surrounding fluid is ac¬ 
tually stable. 

A qualitative explanation of the situation is provided 
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FIG. 9: Schematic plot of the chemical potential p versus 
density p for a system undergoing a liquid-solid transition in 
a finite box volume Vbox with periodic boundary conditions. 
Due to interfacial effects, non-negligible in finite systems, the 
isotherm deviates from p = p CO ex in the two-phase coexis¬ 
tence region, pf < p < p m . The sharp kinks of the curve are 
in reality rounded due to fluctuations as long as Vb ox is finite 
(for Vbox oo, however, all special features of this isotherm 
disappear and only the horizontal straight line p = p coex re¬ 
mains) . At the density pi , the system undergoes a transition 
from homogeneous fluid to a situation of two-phase coexis¬ 
tence (crystalline nucleus (red particles) surrounded by fluid 
(blue particles), see [121 )• Particles in the center of the in¬ 
terfacial region are shown in green, cf. Fig. [3] AT p = p 2 , 
the shape of the crystalline nucleus changes from spherical 
to cylindrical, while at p = p$ a slab configuration (with two 
planar interfaces) sets in. In the density region where the slab 
configuration prevails, p = p CO e X (as well as p = p CO ex) hold 
also for finite volumes. 


in Fig. [9j where the chemical potential is plotted versus 
density: in a finite system, the region of two-phase coex¬ 
istence actually is reduced in comparison with the bulk. 
The homogeneous liquid in fact loses its thermodynamic 
stability not for p = pf but at the higher density p\. 
This enhanced stability of the fluid for densities where in 
the thermodynamic limit the fluid would be metastable is 
due to the fact that the surface free energy cost associated 
with the formation of a nucleus would be too high. The 
same phenomena is well know for the vapor-liquid transi¬ 
tion [53l - i58] , where the transition at p = pi is termed the 
droplet evaporation-condensation transition [53H551 [80] . 
It has been suggested that (55J [501 Pi ~ Pf K P '^oJ 4 i n 
d = 3 dimensions, while the densities P 2 ,P 3 show only 
insignificant dependence on box size. 

A plot of pressure p versus density looks similar to 
Fig. @ when the plot refers to the majority phase, i.e. 
the liquid phase for pf < p < pd = {pp + p m )/2, and the 
crystal for pd < p < p m - Note that the interchange of 
what majority and minority means, at p = pd, happens 
in the regime where the slab morphology prevails, and 
then pe = Pc = Pcoex- On the other hand, in the regions 
where the fluid coexists with compact nuclei (spherical 
droplets, in the vapor to liquid transition, respectively) 
the pressure inside the nucleus is enhanced in compar- 



FIG. 10: (a) Cutout of the cross section through the system, 

for a plane containing the center of mass of the (compact) 
crystalline nucleus for the choice p = 0.527 of the total pack¬ 
ing fraction and N = 8000. Particles belonging to the crystal 
are shown in red; their regular arrangement can be clearly 
recognized. Particles clearly belonging to the fluid are shown 
in blue, while green particles belong to the interfacial region. 
The classification of solid and liquid can be found in [82] • (b) 
Radial density profiles of the order parameters qi and qs, for 
the case where N = 8000,77 = 0.528 and the packing fraction 
in the liquid could be estimated as pi = 0.512 at a pressure 
of pi = 9.74. 


ison with the pressure n the surrounding phase. While 
in situations of phase coexistence the chemical potential 
is homogeneous (as it must be in equilibrium, when the 
coexisting domains may exchange particles freely), the 
pressure distribution is not homogeneous. Even for the 
slab configuration, where the pressure inside both coex¬ 
isting phases is the same, as stated above, a nontrivial 
variation of the pressure tensor components occurs in the 
interfacial regions. 

In this region of packing fractions p, where a stable 
equilibrium between a compact nucleus and surrounding 
fluid is established we record after every 10 • 10 6 Monte 
Carlo steps the configuration of the crystalline nucleus, 
see Fig. |To[ i for an example, and take an average over 
the radial density distribution of the particles belonging 
to this nucleus around its center of mass (Fig. [lO]b). Of 
course, we do not imply that the shape of the crystalline 
nucleus is on average spherical (although this turns out 
to be a rather good approximation). But the informa¬ 
tion shown in Fig. [TOJqb helps us in each configuration of 
the system that is analyzed, to identify the region of the 
fluid that is sufficiently remote from the crystal (where 
< 74,(76 in Fig. |10[ o are already flat, independent of r, and 
have values typical for the liquid phase), where now mea¬ 
surements of the pressure pi and associated packing frac¬ 
tion peijpt) can be performed. Gratifyingly, these data 
(Fig. Ill are fully compatible with the continuation of 
the liquid branch of the isotherm (Fig. [2]) for pressures 
exceeding p CO ex- This result shows that the liquid sur¬ 
rounding the nucleus is in thermal equilibrium, which is 
a basic condition for the analysis presented in the follow¬ 
ing section. 
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FIG. 11: Pressure pe of the liquid surrounding a droplet 
plotted vs. average packing fraction rj of the colloid parti¬ 
cles in the total simulation volume. Three choices of particle 
number N are included, as indicated. For clarity, the region 
of interest is shown also on a plot with strongly magnified 
abscissa scale. Also data for 0.54 < p < 0.575 are included, 
where the crystal forms a slab configuration (cf. Fig. 9), and 
the pressure agrees (apart from fluctuations) with p coe x as it 
should be. 


ESTIMATION OF THE SURFACE EXCESS FREE 
ENERGY OF CRYSTALLINE NUCLEI 
SURROUNDED BY LIQUID FROM A 
THERMODYNAMIC ANALYSIS 


Fig- [Tl] now shows our data for the pressure of the 
fluid phase surrounding the nucleus for the three choices 
of particle number N = 6000, 8000 and 10 000, in the 
regime 0.525 < p < 0.530, where stable nuclei occur. We 
now use a simple extension of the lever rule for phase 
coexistence, adapted to account for finite size effects, to 
infer the volume V* of the crystalline nucleus, namely 

riVbox = T)t(pt)(Vbox ~ V*) + T] c {p c )V* (19) 


Note that 77 is the chosen control variable and when N is 
chosen also Vb ox is fixed. We observe pe as well as rje{pe) 
directly, as discussed above. Since a direct estimation 
of r) c {Pc) in a small fluctuating crystal nucleus suffers 
from relatively large statistical errors (and perhaps also 
systematic errors, if parts of the interfacial region are 
included), we use the appropriate value r] c (p c ) from the 
bulk crystal branch of the isotherm (Fig. [ 2 ]). To find 
Vc{Pc) we use the fact that we know p coex to carry out 
the thermodynamic integration. 

Pc 

Pc{Pc) = Pcoex + ^ J dp/rj c (p) ( 20 ) 

Pcoex 


to find the chemical potential of the crystal branch of 
the isotherm. From Figs. 2|11 we note that r] c (p ) is a 
slightly nonlinear function of p near p CO ex, unlike i)e{p ), 


and therefore it is ne cess ary to go beyond the simple lin¬ 
ear expansion (cf. Eq 15) p c (p c ) ~ /Wx+f ^(Pc~Pcoex) 
(note ri m = rj c {PcoexH The crystal pressure p c then fol¬ 
lows from Eq. [9] since we have determined /1 and pe in 
the liquid region explicitly, the equality of the chemi¬ 
cal potentials determines also p c through Eq. [20] Us¬ 
ing then the data of Fig. [2] where the pressure of the 
crystal p c is determined as function of its packing frac¬ 
tion, we obtain p c (p c ) once p c has been obtained. Then 
the crystallite volume V* is obtained from Eq. [19] with¬ 
out the need of assuming anything about the shape of 
this crystal nucleus. If we would assume that the vol¬ 
ume is strictly spherical, V* = 47 ri?* 3 / 3 , we could try 
to use the radial order parameter profiles of <74 (r) and 
qe(r), Fig. 10a, to estimate R* from the inflection point 
of these curves. Fig. [12] shows effective radii obtained 
in this way, and compares them with an effective radius 
R* e ff = ( 3 V r */ 47 r ) 1 / 3 extracted from the general method 


(using Eqs. 19|20| with the spherical approximation. Un¬ 
expectedly, these data indicate a rather good agreement, 
indicating that for the present model anisotropy effects 
of the interface tension cannot be very important. 

To elaborate the consequences of this finding for the 
theory of homogeneous nucleation of crystals in this sys¬ 
tem, we note that in classical nucleation theory the for¬ 
mation free energy of a nucleus is written in terms of 
volume and surface terms as m 


A F(V*) = -( Pc -pe)V* + F surf (V*) (21) 


In the description of the surface excess free energy 
F sur f(y*) of the crystal nucleus, the dependence of the 
crystal-liquid interface tension 7 (n) on the orientation of 
the interface normal n relative to the crystal axes 
needs to be taken into account. Only for isotropic 7 the 
nucleus is a sphere of radius R and this surface excess 
free energy simply is F sur f = AttR 2 ^ = A iso 7U 2 / 3 with a 
geometric factor A iso = (367T) 1 / 3 . For crystals, however, 
the term Ai SO 7 has to be replaced by a more complicated 
expression 

F SU rf(V) = V 2/3 J 7 (n)ds = A w jV 2/3 , (22) 

A w 


where A w is the surface area of a unit volume whose 
shape is derivable from 7 (n) with the Wulff construction 


and the integration in Eq. 22 is extended over 


this surface. The average interface tension 7 introduced 
in Eq. [22] is defined as 


7 = A, 


-1 


7 (n)ds 


(23) 


The critical nucleus according to the theory of 
homogeneous nucleation satisfies the conditions 
d(AF(V*))/dV* = 0, and hence Eqs. [2l] [22] yield 


V* = 


3 {pc-piY 


11/3 


(24) 
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FIG. 12: (a) Nucleus radius R* extracted from the inflec¬ 

tion point of the order parameter profile qe(r) [within error 
it coincides with inflection point of <74 (r)] plotted vs. average 
packing fraction. 


and the associated free energy barrier is 

A F* = lA w jV* 2 / 3 = l -{p c - Pi )V* 


(25) 


While in the thermodynamic limit the liquid surround¬ 
ing the critical nucleus is metastable, and actually the 
equilibrium condition d(AF(V*)/dV*) = 0 refers to the 
unstable equilibrium on top of a saddle point in con¬ 
figuration space, we here deal with a stable equilibrium 
between the nucleus and the surrounding nucleus in the 
finite box. The equilibrium condition is the same, how¬ 
ever, because in both cases the nucleus can freely ex¬ 
change particles with the surrounding liquid, and the 
chemical potential in the system is constant. Of course, 
the analysis presented cannot be extended to arbitrary 
small box volumes, since working at fixed packing frac¬ 
tion 77 = (n/6)Na 3 /Vb ox the chemical potential in a small 
volume is not strictly constant, but exhibits statistical 
fluctuations which are not independent from the statis¬ 
tical fluctuations of V *, and hence give rise to finite size 
effects. 

However, an attractive feature of Eq.[25]is that actually 
neither a knowledge of A w nor of 7 is needed to predict 
A F*: knowing V* and the pressure difference p c — pe 
yields directly A F*, and hence the need to estimate 7 (n) 
and carry out the Wulff construction is bypassed. 

Fig. 13 shows a plot of A F* vs. V * 2 / 3 , using our sim¬ 
ulation results for the three different choices of N. Grat- 
ifyingly, the data for different N overlap, there do not 
seem to be significant finite size effects present (at least 
not for the regime of V* > 250 explored here). It also is 
interesting that the data do not indicate significant devi¬ 
ations from the straight line behavior, indicating that we 
get a meaningful estimate of the constant A w 7/3 from 
the slope of this straight line. Unlike related work on 



y.2/3 

FIG. 13: Barrier A F* against homogeneous nucleation (in 
units of ksT) plotted vs. V * 2 / 3 , using the data for V* and 
Pc — Pe from the thermodynamic analysis of the equilibrium 
between nucleus and surrounding liquid, for the data for 0.525 
< 77 < 0.530 and N = 6000, 8000 and 10000. The straight line 
is a fit through the data that includes the origin (A F* = 0 
for V* =0), using an independent estimate for the surface 
tension and assuming a spherical nucleus shape (see text). 


vapor to liquid nucleation in simple fluids [56H58] . cur¬ 
vature corrections (of order F * 1 / 3 in this plot) do not 
seem to play any role here, thus there is also no need to 
consider the distinction between the equimolar dividing 
surface (implicit in Eq. 19) and the surface of tension 
which needed to be considered in the vapor to liquid case 

Ell). 


Since we have already seen that the direct estimation of 
nuclei volume from a spherical approximation (Fig. |12| ) is 
in good agreement with the volume determined from the 
thermodynamic analysis, it is of course tempting to try a 
direct spherical approximation of A F* as well, replacing 
A w in Eq. 23 by A iso = ( 367 t) 1 / 3 ) and using instead of 7 


an estimate for the smallest interface tension [ 81 ], which 
occurs for the close-packed 111 crystal surface, namely 
7111 = 1.013(3), yielding A iso 7/3 ~ 1.633. It turns out 
that the resulting straight line is almost indistinguish¬ 
able from the straight line included in Fig. [13] This good 
agreement can also be seen when one considers the pres¬ 
sure difference A p = p c — p coex as function of 77 (Fig. 14) 
and includes there the prediction appropriate for a spher¬ 
ical nucleus, A p = ( 27 /R*)/(^ L — 1 ) for the individual 
choices of TV, using again 7m as an approximation for 7. 
Although the other interface tensions 7 ioo ; 7 no for the 
present model are about 3% larger than 7 m [81], this 
weak anisotropy of the present model does not lead to 
significant deviations from the spherical approximation, 
for the present choice of model parameters. 
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FIG. 14: Pressure difference pt — p CO ex plotted versus the 
average packing fraction p in the simulation box, for particle 
number N = 6000 (topmost set of data), 8000 (middle set of 
data) and 10 000(lowest set), respectively. Data points are 
from the thermodynamic analysis (green symbols represent 
the formula pt — p coex = ^Pf(pf — Pcoex)), curves show the 
prediction pt —p CO ex = c- ( 27 in/H)/(r? m /r?/ - 1 ), using 7 m « 
1.013 (taken from [81]) and c = 1.07 as indicated. 


CONCLUSIONS 


In this paper we have introduced a smooth variant of 
the Asakura-Oosawa model for colloid-polymer mixtures, 
in the limit of small size ratios where the polymers can 
be integrated out in favor of a short range colloid-colloid 
attraction, to study the phase transition from the fluid 
phase to the face-centered cubic crystal. We have paid 
particular attention to phase coexistence in finite vol¬ 
umes, and the information on nucleation barriers that 
one can extract from an analysis of phase coexistence. A 
key ingredient of the discussion has been the observation 
that the chemical potential of the system, considered as a 
function of packing fraction at constant volume, exhibits 
a loop (Fig. [ 9 ]). The enhancement of p relative to p coe x 
in the first descending part of this loop, and the associ¬ 
ated enhancement of the pressure pi in the liquid over 
Pcoex , contains information that can be used to estimate 
the surface excess free energy of the crystalline nucleus. 

In order to be able to study this problem with mean¬ 
ingful accuracy, we have first located phase coexistence 
between bulk phases from the standard “interface veloc¬ 
ity” method, where one prepares a “slab configuration” of 
a crystal coexisting with liquid, separated from it by pla¬ 
nar interfaces, in a simulation box of suitably elongated 
shape. Studying the “velocity” with which the volume of 
the box changes with time (due to interface motion, when 
the crystal shrinks for p < p coex or grows for p > p coe x) , 
one can estimate the coexistence pressure Pcoex- Unfortu¬ 
nately, this method due to large statistical fluctuations 
(Fig. [4| and finite size effects (Fig. [5]) requires a ma¬ 
jor computational effort, and alternative more efficient 
methods to accurately estimate p CO ex would be desirable. 


After this part of our work was completed, a method due 
to Pedersen et al. [ 88 ] came to our attention, which could 
be advantageous for this purpose. 

We have also applied a method where the fluid in an 
elongated simulation box is exposed to one wall with a 
very soft potential, enabling the estimation of the chem¬ 
ical potential p also in a situation where in the bulk the 
packing fraction is too high so that the standard Widom 
particle insertion method would not work (Fig. [6|7| . We 
have shown that by this method p can be estimated 
for the liquid branch of the isotherm, even for pressures 
P ’ Pcoex- 

The crystalline nuclei (coexisting with fluid) have been 
studied both by geometrical identification where the 
crystal-fluid interface is located (Fig. [Top] ) using the 
well-known concept of bond orientational order parame¬ 
ters (£ 4 , (fe, cf. Fig. [ 8 ]) to locally distinguish the phases, 
and by thermodynamic analysis (Figs Tl|l3|14 1. We 
show that for the present model the assumption that 
the anisotropy of 7 (n) can essentially be ignored, and 
therefore the nuclei have effectively a spherical shape, 
works very accurately. We stress, however, that we do 
not expect that this conclusion works in general; if the 
anisotropy 7 (n) is more pronounced, the spherical ap¬ 
proximation must fail, but our thermodynamic method 
which only needs the information on the fluid pressure, 
fluid chemical potential, and nucleus volume, still should 
work. 

Such a more anisotropic case is expected for the present 
model at large z p in Eq.JTJr, when is much smaller than 
rj m . However, an explicit study of this case must be left 
to future work. 

An important feature of our work also is the excel¬ 
lent agreement with classical nucleation theory, A F* = 
A w "/V * 2 / 3 /3 (Fig. 13). This conclusion seems to be at 
variance with both experimental and simulation studies 
on nucleation in colloidal suspensions (see [19] for a re¬ 
view). However, we feel that most of this work has been 
done for conditions where {pe~Pf)/pf is much larger than 
in the present work, so actually the nucleation barrier is 
then rather small; our pressure range 9.54 < pa 3 /ksT < 


10.40 for the pressure of the metastable fluid (Fig. 11) 


translates into packing fractions 0.51 < p < 0.52 of 
the fluid surrounding the crystalline nucleus in our box. 
These numbers are close to pf = 0.494, the relative dis¬ 
tance from the transition where freezing sets in is at most 
6%. In contrast, experiments [3] and simulations for hard 
spheres [831485] typically deal with packing fractions in 
the range 0.53 < p < 0.57 (note that within the available 
accuracy pf in our model agrees with the corresponding 
value for hard spheres). Note that this range of the exper¬ 
iments corresponds to huge values of pt~p C oex (cf. Fig.[2|) , 
and correspondingly the estimate 


R* = 2q /[(j>£ - Pcoex) {Pm/Vf - 1)] 


(26) 
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that classical nucleation theory would imply yields 
droplet radii that do not exceed about two colloid di¬ 
ameters, and barriers that are only of the order of 5 to 
10 kgT. It is clear that for such small nuclei the approxi¬ 
mations of the classical nucleation theory are expected to 
fail. On the other hand, for such dense colloidal suspen¬ 
sions the nucleation rate will be strongly affected by the 
slowing down that sets in when the colloidal glass tran¬ 
sition gana is approached. In fact, Radu and Schilling 
m have pointed out a strong and nontrivial effect due 
to the kinetic prefactor of the nucleation rate. Since the 
slowing down associated with the glass transition is not 
a well understood problem [42], studies attempting to 
test nucleation theory should avoid this region, and fo¬ 
cus on the region close to rjf, as done here. Finally, we 
discuss the issue why the spherical approximation is such 
a good approximation to the actual crystal shape in this 
case. The simple explanation is that the anisotropy pf 
the interfacial tension of the crystal fluid interface in this 
model is very weak, the difference between 7m,7no and 
7100 is a few % only. A similar situation also arises in 
the lattice gas model [ 28 ], at temperatures well above 
the roughening temperature. Recent work |&Ql has stud¬ 
ied the angular dependence of the interface tension in 
this model over a wide range of temperatures, showing 
that the dependence on interface tension on interface ori¬ 
entation amounts to a variation of typically a few % as 
well (this anisotropy vanishes at the critical point of this 
model). The corresponding enhancement of the nucle¬ 
ation barrier was found to be of the order of 10 % as 
well [ 28 ]; only at low temperatures, where facetted crys¬ 
tal shapes occur, large derivations from the spherical 
approximations were found [ 28 ]. Thus, we expect that 
the present findings are generic for all systems where the 
anisotropy of the interfacial tension at melting conditions 
amounts only a few %. 
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